Objectives

Implement the Metropolis algorithm for a one-parameter model, and apply it to generate a posterior parameter sample, for the model “survival”.
Calculate a Markov chain sample of the posterior model parameter distribution. Starting from a first guess jump distribution, construct a simple algorithm to improve the jump distribution for higher efficiency.
Use graphical analysis and convergence tests to asses the quality of the Markov chain samples.

1. Implement an elementary Markov Chain sampler ★

In this exercise, you will write your own code to implement the most elementary MCMC sampler, the Metropolis algorithm for one dimension. MCMC samplers can be used to sample from any distribution, and they are not linked to Bayesian concepts per se. In Bayesian inference, however, they are especially useful to sample from the posterior distribution. You will apply your self-coded sampler to sample from the posterior distribution of the one-parameter model “survival”, given the observed data in the file model_survival.csv.

Before you get started

We want to use the model “Survial” from yesterday’s exercise and the data model_survival.csv to test our sampler. The experiment was started with N=30 individuals (this information is not contained in the data set).

Reading and plotting the data

It’s always a good idea to look at the data first.

Hints

R

data.survival <- read.table("../../data/model_survival.csv", header=TRUE)
barplot(data.survival$deaths, names=data.survival$time,
        xlab="Time interval", ylab="No. of deaths",
        main="Survival data")

Julia

using DataFrames
import CSV
using Plots
survival_data = CSV.read("../../data/model_survival.csv",
                         DataFrame);

bar(survival_data.deaths,
    labels = false,
    xlabel = "Times interval",
    ylabel = "No. of deaths",
    title = "Survival data")

Defining the likelihood, prior, and posterior

In the next step we will need the posterior density (or more accurately: a function that is proportional to the posterior density) for the model “survival” because \[ p(\lambda|y)\propto p(\lambda,y) = p(y|\lambda)p(\lambda) \; . \] This means that we need the likelihood function \(p(y|\lambda)\) and the prior \(p(\lambda)\).

Likelihood

R

The likelihood is given below (copied from Exercise 1):

loglikelihood.survival <- function(y, N, t, par){

    ## Calculate survival probabilities at measurement points
    ## and add zero "at time-point infinity":
    S <- exp(-par*c(0,t))
    S <- c(S,0)

    ## Calcute probabilities to die within observation windows:
    p <- -diff(S)

    ## Add number of survivors at the end of the experiment:
    y <- c(y, N-sum(y))

    ## Calculate log-likelihood of multinomial distribution at point y:
    LL <- dmultinom(y, prob=p, log=TRUE)

    return(LL)
}

Julia

The likelihood can be implemented like this, copied from Exercise 1.

using Distributions

function loglikelihood_survival(N::Int, t::Vector, y::Vector, λ)
    S = exp.(-λ .* t)
    S = [1; S; 0]
    p = -diff(S)
    ## Add number of survivors at the end of the experiment:
    y = [y; N-sum(y)]

    logpdf(Multinomial(N, p), y)
end

Note, that we return -Inf if a negative parameter is given. However, we still need a prior distribution for the model parameter \(\lambda\) that encodes that it must be positive.

Prior

The parameter \(\lambda\) must be positive. Furthermore, we know that values around 0.3 are quite plausible, and that values below 0.1 and above 0.7 are pretty unlikely. Which distribution would you choose for the prior?

Implement a function that returns the logarithm of a probability density of the prior distribution for \(\lambda\).

R

You can use a lognormal prior distribution by understanding and filling-in the block below.

# we use a lognormal distribution, which naturally accounts
# for the fact that lambda is positive

logprior.survival <- function(par){

  # a mean and standard deviation that approximate the prior information given above are:
  mu <- ???    # since the distribution is skewed, choose it higher than 0.3 (the mode)
  sigma <- ??? # a value similar to mu will do

  # calculate the mean and standard deviation in the log-space
  meanlog <- log(mu) - 1/2*log(1+(sigma/mu)^2) #
  sdlog <- sqrt(log(1+(sigma/mu)^2))           #

  # Use these to parameterize the prior
  dlnorm(???)
}

Julia

You can use a lognormal prior distribution by understanding and filling-in the block below.

function logprior_survival(λ)
    # a mean and standard deviation that approximate the prior information given above are:
    m = ?       # the distribution is skewed, choose a mean higher than 0.3 (the mode)
    sd = ?      # a value similar to the mean will do

    # calculate the mean and standard deviation in the log-space
    μ = log(m/sqrt(1+sd^2/m^2))
    σ = sqrt(log(1+sd^2/m^2))

    # Use these to parameterize the prior
    return logpdf(Normal(???), ???)
end

Posterior

Finally, define a function that evaluates the density of the log posterior distribution.

This function should take the data from the file model_survival.scv and a parameter value as inputs.

Check what happens if you pass an “illegal”, i.e. negative parameter value. The function must not crash in this case.

R

log.posterior <- function(par, data) {
    ...
}

Julia

function logposterior_survival(λ, data)
    return ...
end

Code your own Markov Chain Metropolis sampler

Implement your own elementary Markov Chain Metropolis sampler and run it with the posterior defined above for a few thousand iterations. You can find a description of the steps needed in Carlo’s slides.

R

#  set sample size (make it feasible)
sampsize <- 5000

# Create an empty matrix to hold the chain:
chain <- matrix(nrow=sampsize, ncol=2)
colnames(chain) <- c("lambda","logposterior")

# Set start value for your parameter:
par.start <- ???

# compute logosterior value of of your start parameter
survival[1,] <- ???

# Run the MCMC:
for (i in 2:sampsize){

  # Propose new parameter value based on previous one
  # (choose a propsal distribution):
  par.prop <- ???

  # Calculate logposterior of the proposed parameter value:
  logpost.prop <- ???

  # Calculate acceptance probability by using the Metropolis criterion min(1, ... ):
  acc.prob <- ???

  # Store in chain[i,] the new parameter value if accepted, or re-use the old one if rejected:
  if (runif(1) < acc.prob){
    ???
   } else {
    ????
   }
}

Julia

# set your sampsize
sampsize = 5000

# create empty vectors to hold the chain and the log posteriors values
chain = Vector{Float64}(undef, sampsize);
log_posts = Vector{Float64}(undef, sampsize);

# Set start value for your parameter
par_start = ???

# compute logosterior value of of your start parameter
chain[1] = ???
log_posts[1] = ???

# Run the MCMC:
for i in 2:sampsize

  # Propose new parameter value based on previous one (choose a propsal distribution):
  par_prop = ???

  # Calculate logposterior of the proposed parameter value:
  logpost_prop = ???

  # Calculate acceptance probability by using the Metropolis criterion min(1, ... ):
  acc_prob = ???

  # Store in chain[i] the new parameter value if accepted,
  # or re-store old one if rejected:
  if rand() < acc_prob # rand() is the same as rand(Uniform(0,1))
    ???
  else # if not accepted, stay at the current parameter value
    ???
  end
end

Experiment with you sampler: - Plot the resulting chain. Does it look reasonable? - Test the effect of different jump distribution. - Find the value of \(\lambda\) within your sample, for which the posterior is maximal. - Create a histogram or density plot of the posterior sample. Look at the chain plots to decide how many the burn-in sample you need to remove.

Solution

R

Defining the likelihood, prior, and posterior

loglikelihood.survival <- function(y, N, t, par){

    ## Calculate survival probabilities at measurement points
    ## and add zero "at time-point infinity":
    S <- exp(-par*c(0,t))
    S <- c(S,0)

    ## Calcute probabilities to die within observation windows:
    p <- -diff(S)

    ## Add number of survivors at the end of the experiment:
    y <- c(y, N-sum(y))

    ## Calculate log-likelihood of multinomial distribution at point y:
    LL <- dmultinom(y, prob=p, log=TRUE)

    return(LL)
}

Using a log normal distribution as prior

logprior.survival <- function(par){

  # a mean and standard deviation that approximate the prior information given above are:
  mu <- 0.6    # since the distribution is skewed, choose it higher than 0.3 (the mode)
  sigma <- 0.7 # a value similar to mu will do

  # calculate the mean and standard deviation in the log-space
  meanlog <- log(mu) - 1/2*log(1+(sigma/mu)^2) #
  sdlog <- sqrt(log(1+(sigma/mu)^2))           #

  # Use these to parameterize the prior
  dlnorm(par, meanlog=meanlog, sdlog=sdlog, log=TRUE)
}

We also define a function for the log posterior

logposterior.survival <- function(par, data) {
    loglikelihood.survival(data$deaths, N=30, data$time, par) + logprior.survival(par)
}

Code your own Markov Chain Metropolis sampler

# Set sample size (make it feasible)
sampsize <- 5000

# Create an empty matrixto hold the chain:
chain <- matrix(nrow=sampsize, ncol=2)
colnames(chain) <- c("lambda","logposterior")

# Set start value for your parameter:
par.start <- c("lambda"=0.5)

# Compute posterior value of first chain by using your start value
# of the parameter, and the logprior and loglikelihood already defined
chain[1,] <- c(par.start,
               logposterior.survival(par.start, data.survival))

# Run the MCMC:
for (i in 2:sampsize){

  # Propose new parameter value based on previous one (think of propsal distribution):
  par.prop <- chain[i-1,1] + rnorm(1, mean=0, sd=0.01)

  # Calculate logposterior of the proposed parameter value (use likelihood and prior):
  logpost.prop <- logposterior.survival(par.prop, data.survival)

  # Calculate acceptance probability:
  acc.prob <- min( 1, exp(logpost.prop - chain[i-1,2]))

  # Store in chain[i,] the new parameter value if accepted, or re-use the old one if rejected:
  if (runif(1) < acc.prob){
    chain[i,] <- c(par.prop, logpost.prop)
  } else {  # if not accepted, stay at the current parameter value
    chain[i,] <-  chain[i-1,]
  }
}

Plot the resulting chain. Does it look reasonable? Have you selected an appropriate jump distribution?

plot(chain[,"lambda"], ylab=expression(lambda), pch=19, cex=0.25, col='black')

Find the value of \(\lambda\) within your sample, for which the posterior is maximal.

# Find the parameter value corresponding to the maximum posterior density
chain[which.max(chain[,2]),]

##       lambda logposterior 
##     0.191718    -9.751701

Create a density plot of the posterior sample using plot(density()), but remove the burn-in first. Inspect the chain in the figure above: how many samples would you discard as burn-in?

# Create a density plot of the posterior sample using `plot(density())` without burn-in
plot(density(chain[1000:nrow(chain),1]), xlab=expression(lambda), main="")

Julia

Defining the likelihood, prior, and posterior

using Distributions

function loglikelihood_survival(N::Int, t::Vector, y::Vector, λ)
    S = exp.(-λ .* t)
    S = [1; S; 0]
    p = -diff(S)
    ## Add number of survivors at the end of the experiment:
    y = [y; N-sum(y)]

    logpdf(Multinomial(N, p), y)
end

## loglikelihood_survival (generic function with 2 methods)


function logprior_survival(λ)
    # a mean and standard deviation that approximate the prior information given above:
    m = 0.5     # since the distribution is skewed, choose it higher than 0.3 (the mode)
    sd = 0.5    # a value similar to the mean will do

    # calculate the mean and standard deviation in the log-space
    μ = log(m/sqrt(1+sd^2/m^2))
    σ = sqrt(log(1+sd^2/m^2))

    # Use these to parameterize the prior
    return logpdf(LogNormal(μ, σ), λ)
end

## logprior_survival (generic function with 1 method)

Note, we check if a parameter value is possible before calling the likelihood function:

function logposterior_survival(λ, data)
    lp = logprior_survival(λ)
    if !isinf(lp)
        lp += loglikelihood_survival(30, data.time, data.deaths, λ)
    end
    lp
end

## logposterior_survival (generic function with 1 method)

Code your own Markov Chain Metropolis sampler

# set your sampsize
sampsize = 5000

# create empty vectors to hold the chain and the log posteriors values
chain = Vector{Float64}(undef, sampsize);
log_posts = Vector{Float64}(undef, sampsize);

# Set start value for your parameter
par_start = 0.5

# compute logposterior value of of your start parameter
chain[1] = par_start
log_posts[1] = logposterior_survival(par_start, survival_data)

# Run the MCMC:
for i in 2:sampsize
    # Propose new parameter value based on previous one (choose a propsal distribution):
    par_prop = chain[i-1] + rand(Normal(0, 0.01))

    # Calculate logposterior of the proposed parameter value:
    logpost_prop = logposterior_survival(par_prop, survival_data)

    # Calculate acceptance probability by using the Metropolis criterion min(1, ... ):
    acc_prob = min(1, exp(logpost_prop - log_posts[i-1]))

    # Store in chain[i] the new parameter value if accepted, or re-store old one if rejected:
    if rand() < acc_prob # rand() is the same as rand(Uniform(0,1))
        chain[i] = par_prop
        log_posts[i] = logpost_prop
    else # if not accepted, stay at the current parameter value
        chain[i] = chain[i-1]
        log_posts[i] = log_posts[i-1]
    end
end

The so called trace-plot can give us an idea how well the chain has converged:

plot(chain, label=false, xlab="iteration", ylab="λ")

Find the parameter value corresponding to the maximum posterior density

chain[argmax(log_posts)]

## 0.19227282763827366

log_posts[argmax(log_posts)]

## -9.619980909799537

A histogram of the posterior. Note, we removed the first 1000 iterations as burn-in:

histogram(chain[1000:end], xlab="λ", label=false)

2. Sample from the posterior for the monod model ★

The main difference of the “Monod” compared to the “Survival” model is that it has more than one parameter that we would like to infer.

You can either generalize our own Metropolis sampler to higher dimensions or use a sampler from a package.

Read data

We use the data model_monod_stoch.csv:

data.monod <- read.table("../../data/model_monod_stoch.csv", header=T)
plot(data.monod$C, data.monod$r, xlab="C", ylab="r", main='"model_monod_stoch.csv"')

Define likelihood, prior, and posterior

Again, we need to prepare a function that evaluates the log posterior. You can use the templates below.

R

# Logprior for model "monod": lognormal distribution

prior.monod.mean <- 0.5*c(r_max=5, K=3, sigma=0.5)
prior.monod.sd   <- 0.5*prior.monod.mean

logprior.monod <- function(par,mean,sd){

    sdlog    <- ???
    meanlog  <- ???

    return(sum(???)) # we work log-space, hence sum over multiple independent data points
}

# Log-likelihood for model "monod"

loglikelihood.monod <- function( y, C, par){

  # deterministic part:
  y.det <- ???

  # Calculate loglikelihood:
  return( sum(???) )
}

# Log-posterior for model "monod" given data 'data.monod' with layout 'L.monod.obs'

logposterior.monod <- function(par) {
  logpri <- logprior.monod(???)
  if(logpri==-Inf){
    return(-Inf)
  } else {
    return(???) # sum log-prior and log-likelihood
  }
}

Julia

include("../../models/models.jl") # load monod model
using ComponentArrays
using Distributions

# set parameters
par = ComponentVector(r_max = 5, K=3, sigma=0.5)

prior_monod_mean = 0.5 .* par
prior_monod_sd = 0.25 .* par

# Use a lognormal distribution for all model parameters
function logprior_monod(par, m, sd)
    m = ...
    σ = ...
    return sum(...) # sum because we are in the log-space
end


# Log-likelihood for model "monod"
function loglikelihood_monod(par::ComponentVector, data::DataFrame)
    y_det = ...
    return sum(...)
end

# Log-posterior for model "monod"
function logposterior_monod(par::ComponentVector)
    logpri = logprior_monod(...)
    if logpri == -Inf
        return ...
    else
        return ...    # sum log-prior and log-likelihood
    end
end

Create initial chain

First run an initial chain, which most likely has not yet a good mixing.

For the jump distribution, assume a diagonal covariance matrix with reasonable values. The diagonal covariance matrices (with zero-valued off-diagonal entries) correspond to independent proposal distribution.

R

The code below uses the function adaptMCMC::MCMC. If the adaptation is set to FALSE, this corresponds to the basic Metropolis sampler. Of course, feel free to use your own implementation instead!

library(adaptMCMC)


## As start values for the Markov chain we can use the mean of the prior
par.start <- c(r_max=2.5, K=1.5, sigma=0.25)


## sample
monod.chain <- adaptMCMC::MCMC(p  = ...,
                               n = ...,
                               init = ...,
                               adapt = FALSE # for this exercise we do not want to use automatic adaptation
                               )
monod.chain.coda <- adaptMCMC::convert.to.coda(monod.chain) # this is useful for plotting

Plot the chain and look at the rejection rate of this chain to gain information if the standard deviation of the jump distribution should be increased (if rejection frequency is too low) or decreased (if it was too high). What do you think?

monod.chain$acceptance.rate

plot(monod.chain.coda)

pairs(modod.chain$samples)

Julia

The package KissMCMC provides a very basic Metropolis sampler. Of course, feel free to use your own implementation instead!

using KissMCMC: metropolis
using Distributions
using LinearAlgebra: diagm

# define a function that generates a proposal given θ:
Σjump = diagm(ones(3))
sample_proposal(θ) = θ .+ rand(MvNormal(zeros(length(θ)), Σjump))

# run sampler
samples, acceptance_ratio, lp = metropolis(logposterior_monod,
                                           sample_proposal,
                                           par;
                                           niter = 10_000,
                                           nburnin = 0);

## for convenience we convert our samplers in a `MCMCChains.Chains`
using MCMCChains
using StatsPlots
chn = Chains(samples, labels(par))

plot(chn)
corner(chn)

Array(chn)                      # convert it to normal Array

Improve jump distribution

The plot of the previous chain gives us some indication how a more efficient jump distribution could look like.

Try to manually change the variance of the jump distribution. Can we get better chains?
Use the previous chain and estimate its covariance. Does a correlated jump distribution work better?

R

The function adaptMCMC::MCMC takes argument scale to modify the jump distribution:

scale: vector with the variances _or_ covariance matrix of the jump
          distribution.

Julia

Adapt the covariance matrix Σjump of the jump distribution that is used in the function sample_proposal.

Residual Diagnostics

The assumptions about the observation noise should be validated. Often we can do this by looking at the model residuals, that is the difference between observations and model prediction.

Find the parameters in the posterior sample that correspond to the largest posterior value (the so called maximum posterior estimate, MAP).
Run the deterministic “Monod” model with this parameters and analyze the residuals.
- Are they normally distributed?
- Are they independent
- Are they centered around zero?
- Do you see any suspicious structure?

R

The function qqnorm and acf can be helpful.

Julia

The function StatsBase.autocor and StatPlots.qqnorm(x) can be helpful.

Solution

R

Read data

data.monod <- read.table("../../data/model_monod_stoch.csv", header=T)
plot(data.monod$C, data.monod$r, xlab="C", ylab="r", main='"model_monod_stoch.csv"')

Define likelihood, prior, and posterior

source("../../models/models.r") # load the monod model

## Logprior for model "monod": lognormal distribution

prior.monod.mean <- 0.5 * c(r_max=5, K=3, sigma=0.5)
prior.monod.sd   <- 0.5 * prior.monod.mean

logprior.monod <- function(par, mean, sd){
    sdlog <- sqrt(log(1+sd*sd/(mean*mean)))
    meanlog <- log(mean) - sdlog*sdlog/2
    return(sum(dlnorm(par, meanlog=meanlog, sdlog=sdlog, log=TRUE)))
}

## Log-likelihood for model "monod"

loglikelihood.monod <- function(y, C, par){

    ## deterministic part:
    y.det <- model.monod(par, C) # defined in `models.r`

    ## Calculate loglikelihood assuming independence:
    return( sum(dnorm(y, mean=y.det, sd=par['sigma'], log=TRUE )) )
}

## Log-posterior for model "monod"

logposterior.monod <- function(par) {
    lp <- logprior.monod(par, prior.monod.mean, prior.monod.sd)
    if(is.finite(lp)){
        return( lp + loglikelihood.monod(data.monod$r, data.monod$C, par) )
    } else {
        return(-Inf)
    }
}

Create initial chain

library(adaptMCMC)

## As start values for the Markov chain we can use the mean of the prior
par.init <- c(r_max=2.5, K=1.5, sigma=0.25)

## sample
monod.chain <- adaptMCMC::MCMC(p = logposterior.monod,
                               n = 5000,
                               init = par.init,
                               scale = c(1,1,1),
                               adapt = FALSE, # for this exercise we do not
                                              # want to use automatic adaptation
                               showProgressBar = FALSE
                               )

##   generate 5000 samples

monod.chain.coda <- adaptMCMC::convert.to.coda(monod.chain) # this is useful for plotting
summary(monod.chain.coda)

## 
## Iterations = 1:5000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##         Mean     SD Naive SE Time-series SE
## r_max 4.3782 0.3556 0.005029         0.1076
## K     1.6396 0.4078 0.005767         0.1032
## sigma 0.4971 0.1132 0.001601         0.0132
## 
## 2. Quantiles for each variable:
## 
##         2.5%    25%    50%    75% 97.5%
## r_max 3.7319 4.0794 4.4018 4.6310 4.918
## K     0.8047 1.3264 1.6156 1.9390 2.298
## sigma 0.3701 0.4206 0.5033 0.5388 0.688

The acceptance rate is very low:

monod.chain$acceptance.rate

## [1] 0.01

Them poor mixing is also visible from the chain plots:

plot(monod.chain.coda)

Here we see a strong correlation between \(K\) and \(r_{max}\):

pairs(monod.chain$samples)

Improve jump distribution

We use the first chain and estimates it’s covariance matrix which we can use as jump distribution in a second run:

Sigma.jump <- cov(monod.chain$samples) * (1/2)^2

monod.chain2 <- adaptMCMC::MCMC(p = logposterior.monod,
                                n = 5000,
                                init = par.init,
                                scale = Sigma.jump,
                                adapt = FALSE,
                                showProgressBar = FALSE
                                )

##   generate 5000 samples

monod.chain.coda2 <- adaptMCMC::convert.to.coda(monod.chain2)
summary(monod.chain.coda2)

## 
## Iterations = 1:5000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##         Mean      SD Naive SE Time-series SE
## r_max 4.3578 0.32224 0.004557       0.021203
## K     1.6116 0.40042 0.005663       0.025768
## sigma 0.4679 0.08606 0.001217       0.005874
## 
## 2. Quantiles for each variable:
## 
##         2.5%    25%    50%    75%  97.5%
## r_max 3.8035 4.1645 4.3493 4.5602 4.9677
## K     0.9269 1.3349 1.5688 1.8627 2.4678
## sigma 0.3541 0.4164 0.4601 0.5038 0.6187

Everything looks much better now!

monod.chain2$acceptance.rate

## [1] 0.493

plot(monod.chain.coda2)

pairs(monod.chain2$samples)

Residual Diagnostics

First we compute the residuals with the ‘best’ parameters

# get maximum of posterior parameters
par.max <- monod.chain2$samples[which.max(monod.chain2$log.p),]

residuals <- data.monod$r - model.monod(par.max, C=data.monod$C)

Plotting the residuals. We hope not to find any structure, but this model is clearly not perfect:

plot(data.monod$r, residuals, xlab="r")
abline(h=0, col=2)

Quantile-quantile plots are great to see if the residuals are normal distributed (which was out assumption for the likelihood functions). Here we see that we have too heavy tails.

qqnorm(residuals)
qqline(residuals)

We also assumed independence of the observations, i.e. we should not see any auto-correlation. This looks good here:

acf(residuals)

Julia

Read data

using DataFrames
import CSV
using Plots

monod_data = CSV.read("../../data/model_monod_stoch.csv", DataFrame)

## 20×2 DataFrame
##  Row │ C        r
##      │ Float64  Float64
## ─────┼──────────────────
##    1 │     0.5  1.24893
##    2 │     1.0  1.69416
##    3 │     1.5  2.54771
##    4 │     2.0  2.07218
##    5 │     2.5  2.67768
##    6 │     3.0  1.32508
##    7 │     3.5  2.74149
##    8 │     4.0  3.31631
##   ⋮  │    ⋮        ⋮
##   14 │     7.0  3.22343
##   15 │     7.5  3.97951
##   16 │     8.0  3.51488
##   17 │     8.5  3.41961
##   18 │     9.0  3.33122
##   19 │     9.5  3.9776
##   20 │    10.0  4.07524
##           5 rows omitted

scatter(monod_data.C, monod_data.r,
        label=false, xlab="C", ylab="r")

Define likelihood, prior, and posterior

include("../../models/models.jl"); # load the definition of `model_monod`

## model_growth

using ComponentArrays

# set parameters
par = ComponentVector(r_max = 5, K=3, sigma=0.5);

prior_monod_mean = 0.5 .* par;
prior_monod_sd = 0.25 .* par;

# Use a lognormal distribution for all model parameters
function logprior_monod(par, m, sd)
    μ = @. log(m/sqrt(1+sd^2/m^2))
    σ = @. sqrt(log(1+sd^2/m^2))
    return sum(logpdf.(LogNormal.(μ, σ), par)) # sum because we are in the log-space
end

## logprior_monod (generic function with 1 method)


# Log-likelihood for model "monod"
function loglikelihood_monod(par::ComponentVector, data::DataFrame)
    y_det = model_monod(data.C, par)
    return sum(logpdf.(Normal.(y_det, par.sigma), data.r))
end

## loglikelihood_monod (generic function with 2 methods)


# Log-posterior for model "monod"
function logposterior_monod(par::ComponentVector)
    lp = logprior_monod(par, prior_monod_mean, prior_monod_sd)
    if !isinf(lp)
        lp += loglikelihood_monod(par, monod_data)
    end
    lp
end

## logposterior_monod (generic function with 1 method)

Create initial chain

We use the KissMCMC package which provides a very basic metropolis sampler.

using KissMCMC: metropolis
using Distributions
using LinearAlgebra: diagm
using MCMCChains
using StatsPlots

# define a function that generates a proposal given θ:
Σjump = diagm(ones(3))

## 3×3 Matrix{Float64}:
##  1.0  0.0  0.0
##  0.0  1.0  0.0
##  0.0  0.0  1.0

sample_proposal(θ) = θ .+ rand(MvNormal(zeros(length(θ)), Σjump))

## sample_proposal (generic function with 1 method)

# run sampler
samples, acceptance_ratio, lp = metropolis(logposterior_monod,
                                           sample_proposal,
                                           par;
                                           niter = 5_000,
                                           nburnin = 0);


# We convert the sampeles in a `MCMCChains.Chains` object to
# get plotting and summary functions
chn = Chains(samples, labels(par))

## Chains MCMC chain (5000×3×1 Array{Float64, 3}):
## 
## Iterations        = 1:1:5000
## Number of chains  = 1
## Samples per chain = 5000
## parameters        = r_max, K, sigma
## 
## Summary Statistics
##   parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   ess_per_sec
##       Symbol   Float64   Float64   Float64    Float64    Float64   Float64       Missing
## 
##        r_max    4.4118    0.3278    0.0583    26.7810    31.2443    1.0854       missing
##            K    1.5629    0.4372    0.0909    27.0321    52.2027    1.0607       missing
##        sigma    0.4755    0.0635    0.0091    48.7522    84.0187    1.0026       missing
## 
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
## 
##        r_max    3.8982    4.1048    4.3944    4.5986    5.0000
##            K    0.9390    1.2791    1.4229    1.8713    2.6331
##        sigma    0.3655    0.4306    0.4603    0.5000    0.6333

plot(chn)

corner(chn)

Improve jump distribution

The corner plot showed that the parameters \(k\) and \(r_{max}\) are clearly correlated and the chain plots that the jump distribution is too large.

We use the previous sample to estimate the covariance matrix which is then used as new jump distribution.

# take covariance for previous chain
Σjump = cov(Array(chn)) * (1/2)^2

## 3×3 Matrix{Float64}:
##   0.0268606    0.0312574   -0.00118401
##   0.0312574    0.0477896   -0.00051581
##  -0.00118401  -0.00051581   0.0010093

sample_proposal(θ) = θ .+ rand(MvNormal(zeros(length(θ)), Σjump))

## sample_proposal (generic function with 1 method)

samples2, acceptance_ratio2, lp2 = metropolis(logposterior_monod,
                                              sample_proposal,
                                              par;
                                              niter = 5_000,
                                              nburnin = 0);

chn2 = Chains(samples2, labels(par))

## Chains MCMC chain (5000×3×1 Array{Float64, 3}):
## 
## Iterations        = 1:1:5000
## Number of chains  = 1
## Samples per chain = 5000
## parameters        = r_max, K, sigma
## 
## Summary Statistics
##   parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   ess_per_sec
##       Symbol   Float64   Float64   Float64    Float64    Float64   Float64       Missing
## 
##        r_max    4.4069    0.3142    0.0241   197.3471   225.0370    1.0346       missing
##            K    1.6596    0.4260    0.0339   197.0109   150.6084    1.0447       missing
##        sigma    0.4671    0.0657    0.0054   160.3044   207.9493    1.0091       missing
## 
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
## 
##        r_max    3.8331    4.1988    4.3942    4.6027    5.0620
##            K    0.9375    1.3500    1.6178    1.9221    2.5883
##        sigma    0.3585    0.4204    0.4609    0.5077    0.6113

plot(chn2)

corner(chn2)

Much better!

Residual Diagnostics

First we compute the residuals with the ‘best’ parameters

# get maximum of posterior parameters
par_max = samples2[argmax(lp2)]

# compute residuals
residuals = monod_data.r .- model_monod(monod_data.C, par_max)

Plotting the residuals. We hope not to find any structure, but this model is clearly not perfect:

p = scatter(monod_data.r, residuals, label=false, xlab="r", ylab="residual");
hline!(p, [0], label=false)

Quantile-quantile plots are great to see if the residuals are normal distributed (which was out assumption for the likelihood functions). Here we see that we have too heavy tails.

qqnorm(residuals, qqline = :R)               # a bit heavy-tailed

We also assumed independence of the observations, i.e. we should not see any auto-correlation. This looks good here:

using StatsBase
acf = StatsBase.autocor(residuals);
bar(0:(length(acf)-1), acf, label=false, xlab="lag", ylab="correlation")

3. Sample from the posterior for the growth model

The goal is to sample for the posterior distribution of the “Growth” model. This task is very similar to the previous one except that we provide less guidance.

Read data model_growth.csv
Define likelihood, prior, and posterior. What assumptions do you make?
Run MCMC sampler
Check convergence and mixing. If needed, modify the jump distribution
Perform residual diagnostics. Are the assumption of the likelihood function ok?

Exercise 2: Bayesian inference with elementary Markov Chain sampling

Eawag Summer School in Environmental Systems Analysis

1. Implement an elementary Markov Chain sampler ★

Before you get started

Reading and plotting the data

Hints

R

Julia

Defining the likelihood, prior, and posterior

Likelihood

R

Julia

Prior

R

Julia

Posterior

R

Julia

Code your own Markov Chain Metropolis sampler

R

Julia

Solution

R

Defining the likelihood, prior, and posterior

Code your own Markov Chain Metropolis sampler

Julia

Defining the likelihood, prior, and posterior

Code your own Markov Chain Metropolis sampler

2. Sample from the posterior for the monod model ★

Read data

Define likelihood, prior, and posterior

R

Julia

Create initial chain

R

Julia

Improve jump distribution

R

Julia

Residual Diagnostics

R

Julia

Solution

R

Read data

Define likelihood, prior, and posterior

Create initial chain

Improve jump distribution

Residual Diagnostics

Julia

Read data

Define likelihood, prior, and posterior

Create initial chain

Improve jump distribution

Residual Diagnostics

3. Sample from the posterior for the growth model

Solution

R

Read data

Define likelihood, prior, and posterior

Run MCMC sampler

Check convergence and mixing

Residual diagnostics

Julia

Read data

Define likelihood, prior, and posterior

Run MCMC sampler

Check convergence and mixing

Residual diagnostics